function plotrocks(filenumber)
figure(2); clf; hold on;
figure(1); clf; hold on;
surfacehist=[]; timehist=[]; temphist=[]; qhist=[]; avhist=[];
for count=filenumber
    N=readdata(count,1);
    S=readdata(count,2);
    T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
    time=readdata(count,6);
%     dtnom=readdata(count,7);
    dx=readdata(count,8);
    dz=readdata(count,9);
    mx=readdata(count,10);
    mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
    mtype=readdata(count,14);
    mTemp=readdata(count,15);
    %mstressxx=readdata(count,16);
    %mstressxz=readdata(count,17);
    %mstrainxx=readdata(count,18);
    %mstrainxz=readdata(count,19);
    %mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
%     P=readdata(count,23);
Nsurfaces=readdata(count,24);
surfaces=readdata(count,25);

s_focus=round(S/2); depthrange=(25e3:25e3:200e3)+40e3;

count
x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end

[gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1); ind1=find(gridtemp<1);
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); localtemp=interp2(Xgrid',Zgrid',gridtemp',x(s_focus),depthrange);
%temphist=[temphist localtemp];
temphist=[temphist; localtemp'];

mk=zeros(1,N);
mk(find(mtype==1))=3; mk(find(mtype==2))=3; mk(find(mtype>2))=4.08*power((298./(mTemp(find(mtype>2))+273)),0.406);
[gridk,xgrid,zgrid]=bilingrid(mk,mx,mz,dx,dz,1,1);
[gridtemp,xgrid,zgrid]=bilingrid(mTemp,mx,mz,dx,dz,1,1);
tempslice=gridtemp(s_focus,:);
gridq=zeros(S,T);
for s=1:S
    for t=2:T-1
        gridq(s,t)=gridk(s,t)*(gridtemp(s,t+1)-gridtemp(s,t-1))/(dz(t-1)+dz(t));
    end
end
clear('mTemp'); clear ('mtype');
[Zgrid,Xgrid] = meshgrid(zgrid,xgrid); localq=interp2(Xgrid',Zgrid',gridq',x(s_focus),depthrange);
qhist=[qhist; localq'];


sealevel=0.5*(surfaces(1,Nsurfaces)+surfaces(S,Nsurfaces));
surfacehist=[surfacehist surfaces(s_focus,4)];
timehist=[timehist time];
avhist=[avhist trapz(z,tempslice)/z(end)];

startpos=30;
if count>200 & mod(count,50)==0
    figure(2); clf; hold on;
    for i=1:length(depthrange)
        subplot(length(depthrange),1,i);
        plotyy(timehist(startpos:end)/(60*60*365*24),surfacehist(startpos:end)-surfacehist(startpos),timehist/(60*60*365*24),temphist(:,i));
        set(gca,'ydir','rev'); title(['Depth ' num2str(depthrange(i)*1e-3-40) ' km ']);
    end

    figure(3); clf; hold on;
    for i=1:length(depthrange)
        subplot(length(depthrange),1,i);
        plotyy(timehist(startpos:end)/(60*60*365*24),surfacehist(startpos:end)-surfacehist(startpos),timehist/(60*60*365*24),qhist(:,i));
        set(gca,'ydir','rev'); title(['Depth ' num2str(depthrange(i)*1e-3-40) ' km ']);
    end
end

if count>200 & mod(count,50)==0
    figure(1); clf; hold on;
    dtspectral=time/(length(timehist)-1);
    NFFT = 2^nextpow2(length(timehist));
    f=1/dtspectral/2*linspace(0,1,NFFT/2); period=[inf 1./f(2:end)];
    surfaceint=interp1(timehist,surfacehist,0:dtspectral:time); surfaceint(1)=surfacehist(1); surfaceint=surfaceint-surfaceint(1);
    erfinput=1:length(timehist);
    surfaceint=mean(surfaceint)+(surfaceint-mean(surfaceint)).*0.5.*(erf((erfinput-100)/20)-erf((erfinput-(length(timehist)-100))/20));
    surfacespectrum=fft(surfaceint,NFFT)/length(timehist); surfacespectrum=2*abs(surfacespectrum(1:NFFT/2));
    for i=1:length(depthrange)
        subplot(length(depthrange),1,i);
        tempint=interp1(timehist,temphist(:,i),0:dtspectral:time); tempint(1)=temphist(1,i);
        tempint=mean(tempint)+(tempint-mean(tempint)).*0.5.*(erf((erfinput-100)/20)-erf((erfinput-(length(timehist)-100))/20));
        tempspectrum=fft(tempint,NFFT)/length(timehist); tempspectrum=2*abs(tempspectrum(1:NFFT/2));
        %plot(period(2:end)/(60*60*365*24),tempspectrum(2:end),'g'); hold on;
        %plot(period(2:end)/(60*60*365*24),surfacespectrum(2:end));
        [ax,h1,h2]=plotyy(period(2:end)/(60*60*365*24),surfacespectrum(2:end),period(2:end)/(60*60*365*24),tempspectrum(2:end));
        set(ax(1),'xscale','log','yscale','lin'); set(ax(2),'xscale','log','yscale','lin');
        title(['Depth ' num2str(depthrange(i)*1e-3-40) ' km ']);
    end
end

if count>200 & mod(count,50)==0
    figure(4); clf; hold on;
    if count>1
        plotyy(timehist/(60*60*365*24),surfacehist-surfacehist(1),timehist/(60*60*365*24),avhist); set(gca,'ydir','rev');
    end
end

if mod(count,100)==0
    save timeseries depthrange timehist surfacehist temphist qhist avhist;
end
    
% [ax,h1,h2]=plotyy(timehist/(60*60*365*24),surfacehist-surfacehist(1),timehist/(60*60*365*24),temphist(:,1));
% hold(ax(2),'on');
% for i=2:length(depthrange)
%     plot(ax(2),timehist/(60*60*365*24),temphist(:,i));
% end
% axis1=axis(ax(1));
% axis(ax(2),[axis1(1:2) 1000 1600]);

%plotyy(timehist/(60*60*365*24),surfacehist-surfacehist(1),timehist/(60*60*365*24),temphist);
%plot(timehist/(60*60*365*24),surfacehist/surfacehist(1));
%plot(timehist/(60*60*365*24),temphist/temphist(1),'g'); set(gca,'ydir','rev');
pause(0.1);


end
end